library(tidyverse)
library(wesanderson)
library(MetBrewer)
library(sp)
library(sf)
library(classInt)


    #read in world shapefile
    world <- read_rds("figdat/world_borders_WB.rds")
    
    #bring in coefficients collapsed to country level
    coef <- read_csv("figdat/country_coefficients_collapsed.csv")[,c("country","mean_beta","mean_ss_beta")]
    
    #join 
    world$plotOrder <- 1:nrow(world)
    world@data <- left_join(world@data, coef) %>% arrange(plotOrder)
    
    
    
      # Fig 2 panels b and c
      
      pdf("figures/Fig2-bc-raw.pdf", width = 7.2, height = 2*4)
      
              par(mar = c(0,0,0,0))
              par(mfrow = c(2,1))
              
              #panel a - map all
                #define colors
                pal <- met.brewer("Benedictus") 
                int <- classIntervals(world@data$mean_beta, style = "fixed", fixedBreaks = c(-20, seq(-4, 4, 0.5), 20))
                col <- findColours(int, pal)
                
                #plot
                par(bg = 'white')
                plot(world, col = 'gray90', border = NA, ylim = c(-85, 75), xlim = c(-170, 170))
                plot(world, border = NA, col = col, add =T)
                text(x = coordinates(world)[,1], y = coordinates(world)[,2], label = world@data$country, cex = 0.2, col = 'gray90')
              
              
              
              #panel b-  map with SS removed
                  
                #define pal
                pal <- met.brewer("Benedictus") 
                int <- classIntervals(world@data$mean_ss_beta, style = "fixed", fixedBreaks = c(-20, seq(-4, 4, 0.5), 20))
                col <- findColours(int, pal)
                  
                #plot
                par(mar = c(0,0,0,0))
                
                par(bg = 'white')
                plot(world, col = 'gray90', border = NA, 
                     ylim = c(-85, 75), xlim = c(-170, 170))
                plot(world, border = NA, col = col, add =T)
                text(x = coordinates(world)[,1], 
                     y = coordinates(world)[,2], 
                     label = world@data$country, cex = 0.2, col = 'gray90')
                
                
            
      
      dev.off()
      
   








